#include "sgEigen.h"
//#include "Eigen/SuperLUSupport"

using namespace SG::Algebra;
namespace SG::Algebra
{
    template <typename TDataType>
    struct TEigenLLTPrivate
    {
        Eigen::SimplicialLLT<Eigen::SparseMatrix<TDataType>> llt;
    };
}

//
// Containers
//
// The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module;
// it offers high performance and low memory usage. It implements a more versatile variant of
// the widely-used Compressed Column (or Row) Storage scheme.
//
// The case where no empty space is available is a special case, and is referred as the compressed mode.
// It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).
//
// It is worth noting that most of our wrappers to external libraries requires compressed matrices
// as inputs.
//
// The results of Eigen's operations always produces compressed sparse matrices. On the other hand,
// the insertion of a new element into a SparseMatrix converts this later to the uncompressed mode.
//
// see https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
//

// The important point to have in mind when working on sparse matrices is how they are
// stored : i.e either row major or column major. The default is column major.
// Most arithmetic operations on sparse matrices will assert that they have the same
// storage order.
//
// see https://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html
//

//=============================================================================
// TEigenLLT

template <typename TDataType>
TEigenLLT<TDataType>::TEigenLLT ()
{
    d = new TEigenLLTPrivate<TDataType> ();
}

template <typename TDataType>
TEigenLLT<TDataType>::~TEigenLLT ()
{
    if (d) 
    {
        delete d;
    }
}

template <typename TDataType>
void TEigenLLT<TDataType>::compute (const SparseMatrixType& matrix)
{
    d->llt.compute (matrix);
}

template <typename TDataType>
void TEigenLLT<TDataType>::solve (const VectorType& b, VectorType& x) const
{
    x = d->llt.solve (b);
}

template <typename TDataType>
void TEigenLLT<TDataType>::solve (const SparseMatrixType& B, SparseMatrixType& X) const{

    printf("\nERROR Executed a null function 'void TEigenLLT<TDataType>::computeInverse'! \n");
    exit(4);
}

template <typename TDataType>
Info_t TEigenLLT<TDataType>::info () const
{
    // assert (0 && "TEigenLLT<TDataType>::info NOT IMPLEMENTED");
    // return (Info_t) d->llt.info ();

    return (Info_t)0;
}

//=============================================================================
// Explicit instantiations to hide the underlying implementation.
template class ALGEBRA_EXPORT TEigenLLT<Real_t>;
template class ALGEBRA_EXPORT TEigenLLT<Complex_t>;

//=============================================================================
// Explicit specilizations

template <>
const char* ALGEBRA_EXPORT TEigenLLT<double>::getClassName ()
{
    return "TEigenLLT<Real_t>";
}

template <>
const char* ALGEBRA_EXPORT TEigenLLT<Complex_t>::getClassName ()
{
    return "TEigenLLT<Complex_t>";
}